
;;=================================================================
;; Load the GRS thang...
angsep = 10.
;; Get galactic params
defsysv, '!MW', exists = exists
IF ~ exists THEN galactic_params 
d = (n_elements(dvec) NE 0) ? dvec : $
    dindgen(!MW.NBINS)*!MW.BINSIZE + !MW.BINSTART

restore,'./ancillary/grs_kda_data.sav',/ver
restore,'./emaf_paper/bgps_emaf.sav',/ver
n = n_elements(emaf)

restore,'./ancillary/grs_kda_data.sav',/ver
restore,'./emaf_paper/bgps_emaf.sav',/ver
n = n_elements(emaf)

;; Create arrays used for GRS identification
dsep = fltarr(n)
vsep = fltarr(n)
grsi = lonarr(n)

;; Loop through EMAF to check for GRS co-incidence
message,'Checking for GRS co-incidence...',/inf
FOR i=0L, n-1 DO BEGIN
   gcirc, 2, emaf[i].l, emaf[i].b, grs.glon, grs.glat, dis
   mind = min(dis,gind)
   vsep[i] = abs(emaf[i].vlsr - grs[gind].vlsr)
   
   ;; A DSEP =-1 means there is no coincidence for this object
   dsep[i] = (vsep[i] LE grs[gind].delv*1.0) ? mind/60. : -1
   IF dsep[i] GT angsep THEN dsep[i] = -1
   grsi[i] = gind
ENDFOR

emaf.grs = dsep GE 0.
message,'There are '+string(total(emaf.grs),format="(I0)")+$
        ' GRS-associated sources.',/inf

;; Load emaf
emaf.grs_clemdist = grs[grsi].dis
emaf.grs_v        = grs[grsi].vlsr
emaf.grs_l        = grs[grsi].glon
emaf.grs_b        = grs[grsi].glat

;; SELECT ONLY GRS-ASSOCIATED SOURCES
ind = where(emaf.grs, n)
emaf = emaf[ind]

;; Derive GRS distance based on Reid (2009) rotation curve
;; Read in GRS distance correction file
readcol,'./emaf_paper/corrected_grs_distance.txt',bad_gd_cnum,new_gd,$
        format='I,F'
bgdind = WHERE_ARRAY(bad_gd_cnum,emaf.cnum,nbgd)
ngdind = WHERE_ARRAY(emaf.cnum,bad_gd_cnum,nngd)
emaf[bgdind].grs_clemdist = new_gd[ngdind]
grs_dtan = !MW.R0 * cos(emaf.grs_l * !dtor) / cos(emaf.grs_b * !dtor)/1.d3
grs_far = emaf.grs_clemdist GE grs_dtan
print,'N GRS at far KDIST: ',long(total(grs_far))


;; Calculate Reid distances for GRS sources
FOR i=0L,n-1 DO BEGIN
   IF (i+1) MOD 20 EQ 0 THEN $
      print,'Processing object '+string(i+1,format="(I3)")+$
            ' of '+string(n,format="(I0)")
   
   ;; Calculate grs_dist based on GRS KDA resolution
   emaf[i].grs_dist=(grs_far[i]) ? $
                    KDIST(emaf[i].grs_l,emaf[i].grs_b,emaf[i].grs_v,/FAR)/1.d3:$
                    KDIST(emaf[i].grs_l,emaf[i].grs_b,emaf[i].grs_v,/NEAR)/1.d3
ENDFOR

;;=================================================================
;; Back to your regularly-scheduled computations...

nalpha = 251.
alpha = reverse(findgen(nalpha)/(100.)+0.5)
mean_FWHM  = fltarr(nalpha)
mean_ratio = fltarr(nalpha)
mean_pchoo = fltarr(nalpha)
frac_yes   = fltarr(nalpha)
mean_fgy   = dblarr(nalpha)
mean_fgn   = dblarr(nalpha)

grs_dml  = dblarr(nalpha)
grs_dbar = dblarr(nalpha)
dml_c1   = dblarr(nalpha)
dml_c15  = dblarr(nalpha)
dml_c2   = dblarr(nalpha)
dbar_c1  = dblarr(nalpha)
dbar_c15 = dblarr(nalpha)
dbar_c2  = dblarr(nalpha)


ratios = 10.^(findgen(51)/25)
nrat = n_elements(ratios)
pchoo_arr = dblarr(n,nrat)
rhist_arr = dblarr(n,nrat)

indc1  = WHERE( emaf.c GE 0.10,  nc1 )
indc15 = WHERE( emaf.c GE 0.15, nc15 )
indc2  = WHERE( emaf.c GE 0.20,  nc2 )

geq = cgSymbol('geq',/ps)
PREFIX = './emaf_paper/save_files/'

FOR i=0L,nalpha-1 DO BEGIN
   
   message,'Processing for alpha = '+string(alpha[i],format="(F0.2)"),/inf
   
   salpha = string(alpha[i],format="(F0.2)")
   fn = 'irdc_morph_pvec_'+salpha+'.sav'
   IF ~ FILE_TEST(PREFIX+fn) THEN BEGIN
      run_morph_matching,/DRIVEMORPH, ALPHA=alpha[i], /NOPS, /SILENT
      spawn,'mv ./irdc_dist_model/'+fn+' ./emaf_paper/save_files/.'
      restore,PREFIX+fn
      pvec = pvec[where(pvec.has,nhas)]
      save,pvec,filename=PREFIX+fn
   ENDIF
   restore,PREFIX+fn,/ver
   restore,PREFIX+'morph_distances_a'+salpha+'.sav',/ver
   
   
   ;; Set up arrays according to contrast cuts 
   dml = dml[ind,*]
   dbar = dbar[ind,*]
   dpdf = tdpdf[ind,*]
   pvec = pvec[ind]
   
   dml_gy  = where( abs(dml[*,0]  - emaf.grs_dist) LE 1., nml)
   dml_gn  = where( abs(dml[*,0]  - emaf.grs_dist) GT 1.)
   dbarind = where( abs(dbar[*,0] - emaf.grs_dist) LE 1., nbar)
   grs_dml[i]  = double(nml)  / double(n)
   grs_dbar[i] = double(nbar) / double(n)
   
   
   dmlind  = where( abs(dml[indc1,0]  - emaf[indc1].grs_dist) LE 1., nml)
   dbarind = where( abs(dbar[indc1,0] - emaf[indc1].grs_dist) LE 1., nbar)
   dml_c1[i]  = double(nml)  / double(nc1)
   dbar_c1[i] = double(nbar) / double(nc1)
   
   dmlind  = where( abs(dml[indc15,0]  - emaf[indc15].grs_dist) LE 1., nml)
   dbarind = where( abs(dbar[indc15,0] - emaf[indc15].grs_dist) LE 1., nbar)
   dml_c15[i]  = double(nml)  / double(nc15)
   dbar_c15[i] = double(nbar) / double(nc15)
   
   dmlind  = where( abs(dml[indc2,0]  - emaf[indc2].grs_dist) LE 1., nml)
   dbarind = where( abs(dbar[indc2,0] - emaf[indc2].grs_dist) LE 1., nbar)
   dml_c2[i]  = double(nml)  / double(nc2)
   dbar_c2[i] = double(nbar) / double(nc2)
   
   ;;======================================================
   ;; Alpha-testing junk
   nhas = n_elements(pvec)
   FWHM     = fltarr(nhas)
   ratio    = fltarr(nhas)
   pchoo    = fltarr(nhas)
   
   dtan = !MW.R0 * cos(emaf.l * !dtor) / cos(emaf.b * !dtor)
   
   message,'Looping over Objects...',/inf
   ;; Loop over the objects...
   FOR j=0L,nhas-1 DO BEGIN
      ;; Calculate FWHM
      yfit = mpfitpeak(d,pvec[j].invx,A,nterms=4)
      FWHM[j] = 2.355 * A[2] / 1.d3
      
      ;; Calculate Peak Probabilities on either side of d_tan
      totprob = dpdf[j,*]
      totprob /= max(totprob,maxind)
      
      nind = WHERE(d LE dtan[j],nnear)
      find = WHERE(d GT dtan[j],nfar)
      ;;print,nnear,nfar,dtan[j]
      nmax = max(totprob[nind])
      fmax = max(totprob[find])
      
      totprob /= total(totprob)
      
      IF nmax GE fmax THEN BEGIN
         ratio[j] = nmax/fmax 
         pchoo[j] = total(totprob[nind])
      ENDIF ELSE BEGIN
         ratio[j] = fmax/nmax 
         pchoo[j] = total(totprob[find])
      ENDELSE
      
   ENDFOR
   
   mean_FWHM[i]  = median(FWHM)
   mean_ratio[i] = median(ratio)
   mean_pchoo[i] = median(pchoo)
   
   mean_fgy[i] = median(FWHM[dml_gy])
   mean_fgn[i] = median(FWHM[dml_gn])
   
   
;;========================================================
;; Plotting
   
myps,'./emaf_paper/plots/alpha_vs_grs.eps',xsize=13
multiplot,[3,1]

plot,grs_dml,alpha,ytit=cgSymbol('beta'),thick=3,yr=[0.5,3.0],$
     xtit='GRS Matching Fraction for d!dML!n',xr=[0,0.99999],/xst

cgOplot,dml_c1, alpha,thick=3,color='Dodger Blue'
cgOplot,dml_c15,alpha,thick=3,color='Firebrick'
cgOplot,dml_c2, alpha,thick=3,color='Forest Green'

al_legend,/top,/left,charsize=0.8,linestyle=0,box=0,thick=3,$
          color=['Opposite','Dodger Blue','Firebrick','Forest Green'],$
          ['All','C '+geq+' 0.1','C '+geq+' 0.15','C '+geq+' 0.2']

multiplot

plot,grs_dbar,alpha,thick=3,yr=[0.5,3.0],$ ;ytit=cgSymbol_extra('Ypsilon'),$
     xtit='GRS Matching Fraction for d!dbar!n',xr=[0,0.99999],/xst

cgOplot,dbar_c1, alpha,thick=3,color='Dodger Blue'
cgOplot,dbar_c15,alpha,thick=3,color='Firebrick'
cgOplot,dbar_c2, alpha,thick=3,color='Forest Green'

al_legend,/top,/left,charsize=0.8,linestyle=0,box=0,thick=3,$
          color=['Opposite','Dodger Blue','Firebrick','Forest Green'],$
          ['All','C '+geq+' 0.1','C '+geq+' 0.15','C '+geq+' 0.2']


multiplot

;; Figure out the alpha --> FWHM, phcoo issue

cgplot,charsize=1.0,[0,10],[0.5,3.0],/nodata,xtitle='Median FWHM [kpc]',$
       yr=[0.5,3.0],/yst

vline,4,color=cgcolor('Orange'),thick=20

cgOplot,mean_FWHM,alpha,thick=3,color='Opposite'
cgOplot,mean_fgn,alpha,thick=3,color='Firebrick'
cgOplot,mean_fgy,alpha,thick=3,color='Forest Green'

al_legend,/top,/right,box=0,['All','GRS Yes','GRS No'],thick=3,linestyle=0,$
          colors=['Opposite','Forest Green','Firebrick'],charsize=0.85

myps,/done,/mp

   
   
   
   
   
   
   
   
   
   
   
   
   undefine,pvec,dml,dbar
ENDFOR




END
